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Abstract.  The  availability  of  both  global  and  regional  el¬ 
evation  datasets  acquired  by  modern  remote  sensing  tech¬ 
nologies  provides  an  opportunity  to  significantly  improve  the 
accuracy  of  stream  mapping,  especially  in  remote,  hard  to 
reach  regions.  Stream  extraction  from  digital  elevation  mod¬ 
els  (DEMs)  is  based  on  computation  of  flow  accumulation, 
a  summary  parameter  that  poses  performance  and  accuracy 
challenges  when  applied  to  large,  noisy  DEMs  generated 
by  remote  sensing  technologies.  Robust  handling  of  DEM 
depressions  is  essential  for  reliable  extraction  of  connected 
drainage  networks  from  this  type  of  data.  The  least-cost  flow 
routing  method  implemented  in  GRASS  GIS  as  the  mod¬ 
ule  r.  watershed  was  redesigned  to  significantly  improve  its 
speed,  functionality,  and  memory  requirements  and  make  it 
an  efficient  tool  for  stream  mapping  and  watershed  analysis 
from  large  DEMs.  To  evaluate  its  handling  of  large  depres¬ 
sions,  typical  for  remote  sensing  derived  DEMs,  three  differ¬ 
ent  methods  were  compared:  traditional  sink  filling,  impact 
reduction  approach,  and  least-cost  path  search.  The  compari¬ 
son  was  performed  using  the  Shuttle  Radar  Topographic  Mis¬ 
sion  (SRTM)  and  Interferometric  Synthetic  Aperture  Radar 
for  Elevation  (IF S ARE)  datasets  covering  central  Panama  at 
90  m  and  10  m  resolutions,  respectively.  The  accuracy  as¬ 
sessment  was  based  on  ground  control  points  acquired  by 
GPS  and  reference  points  digitized  from  Landsat  imagery 
along  segments  of  selected  Panamanian  rivers.  The  results 
demonstrate  that  the  new  implementation  of  the  least-cost 
path  method  is  significantly  faster  than  the  original  version, 
can  cope  with  massive  datasets,  and  provides  the  most  ac¬ 
curate  results  in  terms  of  stream  locations  validated  against 
reference  points. 
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1  Introduction 

Shuttle  Radar  Topographic  Mission  (SRTM;  Farr  et  al., 
2007)  and  various  airborne  Interferometric  Synthetic  Aper¬ 
ture  Radar  for  Elevation  (IF  S ARE)  surveys  provide  a  new 
generation  of  elevation  data  in  regions  that  have  had  only 
limited,  often  low  resolution  coverage.  These  topographic 
data  sets  are  increasingly  used  to  improve  mapping  of  geo- 
morphic  and  hydrologic  features,  especially  in  remote,  hard 
to  reach  areas  and  at  regional  to  global  scales  (e.g.,  Kinner 
et  al.,  2005;  Lehner  and  Doll,  2004;  World  Wildlife  Fund, 
2009). 

In  spite  of  significant  advances  in  the  development  of  flow 
routing  algorithms  (e.g.  Quinn  et  al.,  1991;  Costa-Cabral  and 
Burges,  1994;  Holmgren,  1994;  Quinn  etal.,  1995;  Tarboton, 
1997),  accurate  extraction  of  drainage  networks  from  Digi¬ 
tal  Elevation  Models  (DEMs)  generated  by  remote  sensing 
technologies  such  as  SRTM  or  IFSARE  remains  challeng¬ 
ing.  To  fully  understand  the  problem,  it  is  important  to  real¬ 
ize  that  the  mapped  elevation  surfaces  include  the  top  of  the 
forest  canopy,  as  well  as  anthropogenic  features,  rather  than 
the  ground  surface  required  for  flow  routing.  Although  ho¬ 
mogeneous  forest  canopy  generally  follows  the  shape  of  the 
ground  surface,  gaps  in  vegetation  that  can  stretch  over  hun¬ 
dreds  of  meters  create  large,  often  nested,  depressions  that 
pose  difficulties  for  flow  routing.  In  addition  to  these  large 
depressions,  the  surface  over  forested  areas  is  noisy,  creating 
numerous  small  depressions  and  barriers  that  further  com¬ 
plicate  drainage  network  extraction.  Therefore,  one  of  the 
important  questions  investigated  in  this  paper  was  whether 
such  data  are  suitable  for  drainage  network  extraction  at  all 
and  if  yes,  how  accurate  are  the  extracted  drainage  networks 
and  which  methods  provide  the  most  reliable  results. 

At  the  same  time,  the  broad  availability  of  elevation  data 
dramatically  increased  the  extent  of  regions  that  can  be  ana¬ 
lyzed  at  relatively  high  resolutions,  given  the  computational 
capabilities  that  can  support  processing  of  massive  DEMs. 
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Therefore,  significant  effort  has  been  devoted  to  the  devel¬ 
opment  of  new  flow  tracing  and  watershed  analysis  algo¬ 
rithms  that  support  efficient  processing  of  large  DEMs  and 
address  the  issue  of  routing  through  complex  nested  depres¬ 
sions  (e.g.  Rivix  Limited  Liability  Company,  2001;  Arge  et 
ah,  2003;  Danner  et  ah,  2007). 

The  most  widespread  method  for  handling  depressions  is 
sink  filling,  up  to  the  level  of  the  sink  spill  point,  com¬ 
bined  with  routing  through  the  resulting  flat  area  (Jenson  and 
Domingue,  1988;  Lairfield  and  Leymarie,  1991;  Planchon 
and  Darboux,  2001;  Wang  and  Liu,  2006).  Improved  sink 
filling  methods  (e.g.  Garbrecht  and  Martz,  1997;  Grimaldi 
et  al.,  2007;  Santini  et  al.,  2009)  first  fill  sinks,  and  then  in¬ 
troduce  a  gradient  to  all  flat  areas  to  provide  non-zero  gra¬ 
dients  for  flow  routing.  The  method  complementary  to  sink 
filling  is  carving  or  breaching  (Rieger,  1998;  Martz  and  Gar¬ 
brecht,  1998)  where  a  channel  is  carved  out  of  each  sink, 
breaking  through  the  (artificial)  obstacle.  Both  principles  can 
be  combined  in  an  impact  reduction  approach  (IRA;  Lindsay 
and  Creed,  2005)  that  for  each  sink  determines  the  method 
that  causes  the  least  impact  on  the  source  dataset.  All  these 
approaches  alter  the  elevation  data  in  order  to  ensure  full 
drainage,  assuming  that  sinks  are  artifacts  created  either  by 
too  low  (artificial  pits)  or  too  high  (artificial  drainage  blocks) 
elevation  values. 

An  alternative  to  the  modification  of  elevation  data  is  to 
determine  the  least-cost  drainage  paths  (LCP)  through  un¬ 
altered  terrain  and  out  of  sinks.  LCP  search  methods  were 
generally  designed  to  find  the  shortest  or  fastest  route  from  a 
starting  point  to  a  given  destination,  used  for  example  in  car 
navigation  devices  or  to  generate  cumulative  cost  surfaces.  A 
particular  LCP  search  method  (A*  Search;  Hart  et  al.,  1968) 
was  adapted  for  flow  routing  and  watershed  analysis  and  im¬ 
plemented  as  the  module  r.watershed  (Ehlschlaeger,  1989)  in 
GRASS  GIS  (Neteler  and  Mitasova,  2008,  see  GRASS  De¬ 
velopment  Team  2010  for  binaries  and  source  code).  When 
applied  to  ILSARE  derived  DEMs,  this  approach  has  pro¬ 
vided  more  accurate  flow  routing  through  large,  nested  de¬ 
pressions  with  fewer  artifacts  than  traditional  sink  filling 
(Kinner  et  al.,  2005).  However,  the  module  was  not  de¬ 
signed  to  handle  the  massive  DEMs  that  are  currently  avail¬ 
able  (Arge  et  al.,  2003).  To  render  the  module  suitable  for 
analysis  of  very  large  datasets  at  ever  higher  resolutions,  and 
to  reduce  the  artifacts  in  flow  patterns  due  to  the  single  flow 
direction  (SLD)  algorithm  (Quinn  et  al.,  1991;  Holmgren, 
1994;  Tarboton,  1997),  redesign  of  the  implementation  and 
addition  of  more  flow  routing  options  was  necessary. 

In  this  study,  we  present  substantial  improvements  to  the 
LCP  method  implementation  and  evaluation  of  its  efficiency 
and  accuracy  when  applied  to  radar-based  DEMs.  We  com¬ 
pare  the  performance  of  the  improved  module  with  (1)  its 
original  implementation,  (2)  a  module  based  on  a  disk  in¬ 
put/output  (I/O)  efficient  method  (Arge  et  al.,  2003),  and 
(3)  the  IRA  implementation  in  the  Terrain  Analysis  System 
(TAS,  Lindsay  and  Creed,  2005).  We  also  analyze  the  impact 


of  mapping  technology  (ILSARE,  SRTM),  resolution,  and 
DEM  resampling  on  the  accuracy  of  the  extracted  drainage 
networks.  Although  the  LCP  method  has  been  used  since 
early  90ies,  the  literature  that  describes  the  algorithm,  its 
implementation  and  properties  is  very  limited  (Ehlschlaeger, 
1989).  This  paper  aims  to  fill  this  gap  in  literature  and  to 
highlight  the  value  of  this  method  and  its  improvements  for 
analysis  of  modem  DEMs. 


2  Methods 

2.1  Improved  least  cost  path  search  algorithm 

The  LCP  algorithm  (A*  Search,  Hart  et  al.,  1968; 
Ehlschlaeger,  1989)  used  for  flow  routing,  flow  accumula¬ 
tion,  and  watershed  analysis  of  raster-based  DEMs  was  re¬ 
designed  to  increase  the  processing  speed  and  decrease  mem¬ 
ory  consumption.  The  implementation  in  the  GRASS  mod¬ 
ule  r.watershed  (Fig.  1)  starts  with  potential  outlet  points. 
Natural  ultimate  outlets  are  e.g.  river  mouths  opening  into 
oceans  or  lakes  without  outflow.  On  gridded  elevation  mod¬ 
els,  potential  outlets  are  grid  cells  along  the  map  boundaries 
or  cells  with  at  least  one  neighbor  with  unknown  elevation, 
e.g.  masked  ocean.  All  potential  outlet  grid  points  are  in¬ 
serted  into  a  list  sorted  by  costs.  Costs  are  measured  as  eleva¬ 
tion  and  order  of  addition  to  the  sorted  list  (grid  cells  added 
earlier  have  higher  precedence  in  case  of  equal  elevation). 
For  the  actual  search,  the  grid  cell  with  the  lowest  eleva¬ 
tion  (smallest  cost)  is  extracted  and  removed  from  the  sorted 
list,  marked  as  processed  and  its  neighbors  are  investigated. 
This  causes  the  search  to  proceed  along  the  least  steep  up¬ 
hill  slope.  At  each  step  during  the  search,  only  neighboring 
grid  cells  that  are  not  yet  in  the  search  list  and  not  yet  pro¬ 
cessed  are  added  to  the  list,  and  drainage  direction  for  these 
neighboring  grid  cells  is  set  towards  the  current  grid  cell.  If 
a  depression  is  encountered,  the  search  follows  the  steepest 
downhill  slope  to  the  bottom  of  a  depression  and  then  pro¬ 
ceeds  again  along  the  least  steep  uphill  slope.  The  search 
proceeds  in  this  manner  upstream  along  the  least  steep  slope 
and  terminates  when  all  grid  points  have  been  processed.  The 
LCP  search  provides  two  results:  (a)  flow  direction  for  each 
cell  in  a  standard  D8  manner  (O’Callaghan  and  Mark,  1984), 
and  (b)  the  order  in  which  cells  must  be  processed  for  flow 
accumulation. 

The  search  component  of  the  original  implementation  was 
modified  by  augmenting  it  with  a  minimum  heap  used  as  pri¬ 
ority  queue  (Atkinson  et  al.,  1986;  Metz  and  Ehlschlaeger, 
2010),  which  led  to  a  substantial  speed  improvement.  In 
the  original  implementation,  the  time  needed  to  keep  the 
sorted  list  sorted  increased  exponentially  with  the  number 
of  grid  cells,  while  in  the  new  implementation  this  time  in¬ 
crease  is  logarithmic.  Removal  of  redundant  information  in 
intermediate  data  as  well  as  reduced  memory  requirements 
for  both  the  search  and  the  flow  accumulation  components 
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Fig.  1.  Flowchart  illustrating  the  least  cost  path  search  algorithm. 


Fig.  2.  SRTM  and  IFSARE  coverage  of  Panama.  GPS  ground  control  points  are  displayed  in  black. 


have  further  enhanced  the  module’s  capacity  to  process  larger 
DEMs.  For  massive  datasets  that  can  not  be  processed  with 
the  amount  of  memory  available,  the  module  can  optionally 
use  external  memory  with  intermediate  data  stored  on  disk. 
These  changes,  aimed  at  improving  the  computational  per¬ 
formance,  did  not  affect  the  quality  of  the  results. 

To  improve  the  flow  routing  accuracy  and  to  reduce  arti¬ 
facts  in  the  flow  accumulation  pattern,  a  multiple  flow  di¬ 
rection  (MFD)  algorithm,  based  on  Holmgren  (1994),  was 
implemented  using  the  order  determined  by  the  LCP  search 
with  an  option  to  control  the  strength  of  flow  convergence. 
For  grid-based  elevation  models,  an  MFD  approach  con¬ 
forms  better  to  theoretical  surface  flow  dispersion  than  single 
flow  direction  in  which  only  one  neighboring  cell  can  receive 
surface  flow  from  the  current  cell  (Freeman,  1991;  Quinn  et 
al.,  1991;  Holmgren,  1994;  Tarboton,  1997).  Initial  D8  flow 
directions  as  determined  by  the  LCP  search  are  used  to  route 
flow  out  of  depressions  when  all  unprocessed  neighboring 
grid  cells  have  higher  elevation  than  the  current  grid  cell. 


2.2  Elevation  data 

Two  different  sources  of  elevation  data  were  used  to  demon¬ 
strate  the  improvements  in  the  LCP  method  (Fig.  2).  Coun¬ 
trywide  elevation  coverage  of  Panama  was  available  as  a 
3  arc  second  resolution  SRTM  DEM  (Farr  et  al.,  2007). 
SRTM  version  2.1  provided  by  the  NASA  Jet  Propul¬ 
sion  Laboratory  (http://www2.jpl.nasa.gov/srtm/,  download 
at  http://dds.cr.usgs.gov/srtm/)  was  selected  for  this  study  as 
the  most  reliable  version  in  terms  of  accuracy  and  minimal 
artifacts,  after  evaluating  properties  of  the  SRTM  products 
available  at  the  time  of  writing  (vl,  v2.1,  v3,  v4.1).  The 
published  absolute  vertical  error  of  the  version  2.1  SRTM 
DEM  is  6.2  m  for  South  America  including  Panama  (Farr 
et  al.,  2007),  although  several  studies  report  higher  accura¬ 
cies  (see  e.g.  Rodriguez  et  al.,  2006).  A  recent  airborne  IF¬ 
SARE  survey  has  provided  new,  more  detailed  information 
about  the  topography  in  central  Panama  at  10  m  resolution, 
using  technology  with  published  vertical  accuracy  of  ~3  m 
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Fig.  3.  Shaded  relief  of  a  Landsat  false  color  composite  (R,G,B  = 
4,5,3)  overlaid  with  GPS  reference  points  in  red  and  points  digitized 
from  Landsat  in  grey. 

(Andersen  et  ah,  2006).  IFSARE  (InterFerometric  Synthetic 
Aperture  Radar  -  Elevation)  is  derived  in  a  way  very  simi¬ 
lar  to  SRTM,  with  the  main  difference  that  IFSARE  is  air¬ 
borne  and  SRTM  data  was  recorded  with  sensors  onboard 
the  Space  Shuttle  Endeavour.  Both  IFSARE  and  SRTM  used 
two  antennas  to  collect  radar  data;  elevation  information  was 
derived  by  analyzing  the  data  recorded  by  the  two  antennas 
(the  general  concept  of  the  new  satellite -borne  TanDEM-X 
mission  is  similar,  collecting  radar  data  with  two  antennas). 
Neither  of  the  surveys  penetrated  the  forest  canopy  to  the 
ground  surface  and  the  elevation  surface  was  over  large  re¬ 
gions  defined  by  a  triple  canopy  tropical  forest  environment 
with  tree  heights  of  more  than  30  m  above  ground.  Given 
the  widespread  use  of  the  SRTM  data,  one  of  the  important 
questions  explored  in  this  paper  was  whether  such  data  are 
suitable  for  stream  extraction  at  all,  and  if  yes,  what  is  the  ac¬ 
curacy  and  which  methods  provide  the  most  reliable  results 
in  this  challenging  environment. 

SRTM  v2  tiles  covering  all  of  Panama  were  combined  and 
gaps  in  the  dataset  filled  using  the  regularized  spline  with  ten¬ 
sion  (RST)  interpolation  method  (GRASS  module  r.fillnulls\ 
Mitasova  and  Mitas,  1993).  The  seamless  SRTM  coverage 
was  then  reprojected  with  bicubic  interpolation  from  geo¬ 
graphic  to  the  UTM  zone  1 7N  coordinate  system  with  90  m 
resolution  to  keep  reprojection  modifications  to  a  minimum. 
To  evaluate  the  impact  of  resampling  on  stream  extraction 
accuracy  (see  e.g.  Valeriano  et  al.,  2006)  and  to  generate 
data  for  additional  performance  testing,  the  reprojected  90  m 
SRTM  DEM  was  reinterpolated  to  30  m  resolution  using  the 
RST  method. 


The  IFSARE  data  were  provided  as  a  10  m  resolution 
DEM  (Kinner  et  al.,  2005)  in  the  UTM  zone  17N  coordi¬ 
nate  system  and  did  not  require  additional  processing.  For 
additional  testing  purposes,  the  IFSARE  DEM  was  down- 
sampled  to  30  m  resolution  using  bicubic  interpolation  in  or¬ 
der  to  have  DEMs  from  two  sources  (SRTM  and  IFSARE) 
with  different  level  of  detail,  but  identical  resolution  of  30  m. 

2.3  Stream  location  data 

Two  sets  of  data  points  were  used  to  evaluate  horizontal  accu¬ 
racy  of  the  drainage  network  extraction  methods:  (i)  stream 
segments  digitized  from  Landsat  imagery  (TM  5,  year  2000, 
scene  id  LT50120542000087XXX02,  provided  by  the  United 
States  Geological  Survey  -  USGS)  and  (ii)  GPS  field  mea¬ 
surements.  Bands  3,  4,  and  5  of  the  Landsat  TM  scene  were 
used  at  its  native  resolution  of  30  m.  The  geographic  projec¬ 
tion  for  this  study  (UTM  17  North)  was  identical  to  the  pro¬ 
jection  of  this  Landsat  scene,  i.e.  the  Landsat  grid  geometry 
and  values  were  not  modified.  A  Landsat  false  color  com¬ 
posite  with  Red,  Green,  and  Blue  assigned  to  the  channels 
4,  5,  and  3  respectively,  clearly  separated  vegetation  from 
waterbodies  and  provided  the  background  for  manual  digiti¬ 
zation  of  points  along  the  river  centerlines.  Only  streams  and 
rivers  with  a  width  less  or  equal  to  4  grid  cells  (120  m)  were 
digitized  in  order  to  exclude  lakes,  anabranching  or  braided 
rivers,  and  rivers  too  broad  to  reliably  determine  a  stream 
center  using  Landsat  imagery  (Fig.  3).  The  horizontal  ac¬ 
curacy  of  these  manually  digitized  reference  points  is  given 
by  the  imagery  resolution  of  30  m  or  better,  depending  on 
the  stream  width  and  the  location  of  the  stream  centerline  in 
relation  to  the  grid  cells. 

GPS  points  were  collected  in  the  field  using  a  Corval¬ 
lis  Microtechnology  March  v3.7  GPS  unit  with  ~2m  posi¬ 
tional  accuracy.  Sites  in  the  Chagres  river  watershed  were 
measured  during  the  years  2002  to  2007  and  sites  at  lower 
reaches  of  most  major  rivers  across  Panama  were  measured 
during  2005  and  2009  years.  GPS  points  were  collected 
along  clearly  identifiable  perennial  rivers  at  locations  where 
a  reliable  GPS  signal  was  available.  The  GPS  measurements 
included  a  larger  proportion  of  points  acquired  in  mountain¬ 
ous  regions  along  smaller  rivers,  while  most  of  the  points 
located  along  larger  rivers  in  low-gradient  flood  plain  and 
coastal  plain  landscapes  were  digitized  from  Landsat  im¬ 
agery  (Fig.  3). 

The  SRTM  water  body  database  was  not  used  for  the  ac¬ 
curacy  evaluation  because  this  cannot  be  regarded  as  an  in¬ 
dependent  data  source  since  it  was  used  in  the  creation  of  the 
SRTM  DEM  v2,  and  because  even  the  larger  rivers  visible  on 
Landsat  imagery  were  still  too  narrow  to  be  captured  in  the 
SRTM  water  body  database. 

We  have  also  considered  and  rejected  the  use  of  hydrog¬ 
raphy  (blue  lines)  from  topographic  maps  because  of  their 
low  reliability  and  a  relatively  coarse  scale  of  1 : 50  000.  Sev¬ 
eral  studies  that  have  compared  the  blue  lines  with  on-ground 
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Fig.  4.  Extraction  of  a  single  main  channel:  (A)  Location  of  areas  shown  in  images  (B),  (C),  and  (D);  (B)  The  starting  points  of  the  extracted 
streams  (red  lines)  were  well  upstream  of  the  reference  points  (black  crosses).  Flow  accumulation  derived  from  IFSARE  10  m  is  used  as 
background.  The  extracted  streams  (red  lines)  are  located  within  broader  stream  tubes  in  low-gradient  areas  close  to  the  coast,  for  both  (C) 
LCP-derived  streams  and  (D)  streams  derived  with  sink  filling.  IFSARE  10  m  was  used  for  both  (C)  and  (D)  representing  the  same  area. 


surveys  of  existing  streams  demonstrated  that  the  blue  lines 
were  the  least  accurate  source  of  stream  location  informa¬ 
tion  (see  e.g.  Colson,  2006  and  the  references  therein),  of¬ 
ten  based  on  old  data  digitized  from  aerial  photographs  with 
large  errors  in  locations  where  the  streams  were  not  visible, 
due  to  the  presence  of  vegetation.  In  fact,  one  of  the  major 
motivations  for  deriving  the  drainage  network  from  the  IF¬ 
SARE  data  in  Panama  was  the  inadequate  scale  of  available 
blue  lines  with  many  larger  streams  missing  and  with  errors 
in  the  mapped  streams  location.  Until  recently,  the  relative 
inaccessibility  of  some  regions  (e.g.  the  upper  Chagres  river) 
and  an  incessant  cloud  cover  have  discouraged  mapping  in 
these  areas  (Kinner  et  al.,  2005). 

2.4  Flow  accumulation  and  evaluation  of  extracted 
drainage  networks 

Several  tests  were  performed  to  evaluate  the  tested  flow  rout¬ 
ing  algorithms.  The  first  test  evaluated  the  processing  speed 
and  capability  to  handle  large  DEMs  when  computing  flow 
accumulation  maps.  In  this  test,  the  performance  of  the  new 
LCP  implementation  was  compared  with  the  old  r.  watershed 
module  (GRASS  Development  Team,  2008)  and  with  a  rel¬ 
atively  new,  I/O-efificient  algorithm,  used  in  the  module 


rterraflow  (Arge  et  al.,  2003).  All  computations  were  done 
on  the  same  hardware  (AMD  Athlon  X2  3  GHz  and  8  GB 
RAM)  but  different  operating  systems  (Linux  64  bit  for  the 
GRASS  modules  r. watershed  and  rterraflow ,  and  Windows 
XP  32  bit  for  TAS  GIS  which  is  a  32  bit  application  for  Mi¬ 
crosoft  Windows  and  was  not  available  as  64  bit  application 
at  the  time  of  writing). 

The  second  test  evaluated  the  horizontal  accuracy  of  the 
extracted  drainage  paths  for  algorithms  with  different  han¬ 
dling  of  depressions.  Three  approaches  were  tested:  (i)  sink 
filling  as  implemented  in  rterraflow,  (ii)  the  impact  reduction 
approach  (IRA)  as  implemented  in  TAS  GIS  (Lindsay  and 
Creed,  2005),  and  (iii)  the  LCP  algorithm  as  implemented  in 
the  enhanced  version  of  r. watershed  in  GRASS  GIS  6.4.  To 
test  the  sink  filling  method,  rterraflow  was  used  because  of 
its  capability  to  efficiently  process  massive  DEMs. 

The  third  test  evaluated  the  impact  of  DEM  resolution 
and  level  of  detail  on  the  horizontal  accuracy  of  the  ex¬ 
tracted  drainage  paths  for  all  tested  algorithms  (; r.watershed ' 
rterraflow ,  and  TAS). 

The  tests  were  performed  for  the  central  Panama  region 
using  both  the  IFSARE  and  SRTM  data  at  the  10  m,  30  m 
and  90  m  resolutions  by  computing  flow  accumulation  with 
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MFD  dispersal  and  extracting  drainage  networks  from  the 
computed  flow  accumulation. 

Drainage  networks  were  extracted  from  all  MFD  flow  ac¬ 
cumulation  maps  using  a  minimum  upstream  catchment  area 
of  100000  m2  as  a  threshold.  This  threshold  was  selected 
for  testing  purposes  and  to  ensure  that  all  measured  or  dig¬ 
itized  points  were  located  downstream  from  the  extraction 
starting  point  and  that  all  points  had  a  drainage  path  associ¬ 
ated  with  it  (Fig.  4b).  Accurate  identification  of  the  channel 
head  zones  was  beyond  the  scope  of  this  paper  because  it  re¬ 
quires  more  than  elevation  data  as  it  is  often  significantly  in¬ 
fluenced  by  local  geology,  groundwater  level,  and  land  cover 
(North  Carolina  Division  of  Water  Quality  NCWQ,  2010). 
According  to  the  NCWQ  (2010)  study,  “the  stream  origins 
usually  occur  as  transition  zones  in  which  the  location  and 
length  of  the  zone  is  subject  to  fluctuations  in  groundwater 
levels  and  precipitation”.  The  tested  90  m  and  30  m  resolu¬ 
tion  SRTM  DEMs  measured  over  at  least  30  m  high  canopy 
of  dense  tropical  forest  did  not  provide  sufficient  informa¬ 
tion  for  accurate  identification  and  mapping  of  the  dynamic 
channel  heads  zones.  Therefore  our  focus  was  on  mapping 
streams  and  rivers  that  have  sufficiently  large  contributing 
areas  that  the  existence  of  a  well  defined  stream  can  be  ex¬ 
pected. 

In  order  to  extract  a  single  main  channel  in  broader  flood- 
plains  where  MFD  flow  accumulation  generates  several  cells 
wide  stream  tubes  that  may  include  cells  above  the  given 
threshold,  a  new  stream  starting  point  was  defined  only  if 
it  was  not  located  within  such  a  stream  tube  (Fig.  4).  Down¬ 
stream  channel  tracing  followed  the  steepest  slope  and  max¬ 
imum  flow  accumulation  within  the  stream  tube.  To  egress 
from  depressions  using  the  LCP  method,  streams  again  fol¬ 
lowed  the  maximum  flow  accumulation  which  was  in  these 
cases  identical  to  the  drainage  direction  initially  determined 
during  the  LCP  search.  As  a  result,  the  drainage  paths  ex¬ 
tracted  from  MFD  flow  accumulation  were  a  single  grid  cell 
wide  and  were  easily  vectorized  for  further  analysis.  This 
approach  of  extracting  single  cell  wide  dendritic  and  topo¬ 
logically  correct  drainage  networks  from  flow  accumulation 
obtained  with  any  kind  of  flow  distribution  method  avoids 
the  problems  associated  with  skeletonizing  rasterized  MFD 
channels.  The  MFD  capabilities  to  map  flow  dispersal  in 
floodplains  can  then  be  preserved  and  single  flow  direction 
(SFD)  routing  is  not  necessary  for  flow  accumulation  and 
subsequent  channel  extraction  (although  it  remains  an  op¬ 
tion). 

Assessment  of  drainage  network  extraction  accuracy 
was  then  performed  by  calculating  distances  between  the 
drainage  paths  derived  from  the  IF  S ARE  and  SRTM  DEMs 
and  the  points  measured  along  the  streams,  consisting  of  338 
on-ground  GPS  measurements  and  995  points  digitized  along 
rivers  from  the  Landsat  imagery  (Figs.  3,4). 

In  order  to  determine  which  method  or  level  of  detail  pro¬ 
vided  more  accurate  results,  we  determined  the  distance  of 
reference  points  to  the  nearest  extracted  drainage  path.  The 


difference  in  the  distance  to  reference  points  between  meth¬ 
ods  or  level  of  detail  was  then  used  to  assess  which  method 
or  level  of  detail  provided  more  accurate  results.  Distance 
differences  were  tested  with  two-tailed  Student’s  t-tests  for 
statistical  significance  and  a  =  0.05.  If  distance  differences 
were  significantly  different  from  zero,  i.e.  one  of  the  tested 
methods/levels  of  detail  was  more  accurate  than  the  other,  the 
sign  of  the  difference  indicated  which  method/level  of  detail 
was  more  accurate. 


3  Results 

3.1  Performance  comparison 

To  demonstrate  the  improvement  in  computational  perfor¬ 
mance  of  the  new  LCP  method  implementation,  the  process¬ 
ing  times  needed  by  the  old  and  new  version  of  r.  watershed , 
and  r.terraflow  were  compared,  with  the  results  summarized 
in  Table  1.  For  flow  accumulation  computation,  the  new 
version  was  350  times  faster  than  the  old  version  for  the 
central  Panama  area  represented  by  a  relatively  small  DEM 
with  2  million  grid  cells  at  90  m  resolution.  The  improve¬ 
ment  was  even  more  dramatic  for  the  countrywide  coverage 
at  90  m  resolution  with  27  million  grid  cells:  the  new  version 
was  about  1750  times  faster.  Although  absolute  processing 
times  are  dependent  on  the  specific  hardware  used  for  test¬ 
ing,  we  expect  that  relative  time  differences  will  be  simi¬ 
lar  on  other  systems  because  the  software  optimizations  are 
hardware-independent.  Processing  the  larger  regions  used 
in  this  study  with  the  old  version  was  impractical  because 
it  could  easily  take  days  (Arge  et  al.,  2003),  whereas  the  new 
version  required  only  minutes  to  execute.  The  computational 
time  needed  by  the  I/O-efihcient  r.terraflow  was  significantly 
lower  than  that  for  the  old  r. watershed ,  but  the  r.terraflow 
module  was  not  as  fast  as  the  new  LCP  implementation  on 
our  test  system  (Table  1).  We  have  also  considered  the  case 
when  the  data  do  not  fit  into  memory  and  r.watershed  uses 
segmented  processing,  with  intermediate  data  being  stored 
on  disk.  This  leads  to  longer  processing  times  than  for  the 
all-in-memory  mode.  Processing  times  observed  on  our  test¬ 
ing  system  in  the  segmented  mode  were  still  shorter  than  for 
r.terraflow ,  which  uses  I/O-efficient  algorithms  specifically 
designed  for  large  datasets  (Arge  et  al.,  2003).  The  size  of 
intermediate  data  created  by  r.watershed  in  the  segmented 
mode  was  less  than  20%  of  those  created  by  r.terraflow.  This 
can  be  of  advantage  if  intermediate  data  are  just  a  bit  too  large 
to  fit  into  memory.  It  is  important  to  note  that,  theoretically, 
r.terraflow  with  its  advanced  I/O-eflicient  algorithms  should 
be  faster  than  r.watershed  for  DEMs  much  larger  than  the 
ones  processed  here  and  on  systems  with  less  physical  mem¬ 
ory. 

TAS  GIS  was  not  included  in  the  performance  compar¬ 
isons,  because  its  memory  requirements  for  processing  the 
10  m  IFSARE  DEM  and  the  30  m  SRTM  DEM  for  all  of 
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Table  1.  Processing  time  for  the  different  DEMs  on  a  Linux  64  bit  system  with  an  AMD  Athlon  X2  3  GHz  CPU  and  8  GB  RAM. 


SRTM  90  m 

Central  Panama 
1.9  million  cells 

IFSARE  30  m 
central  Panama 
17  million  cells 

SRTM  90  m 
all  of  Panama 

27  million  cells 

IFSARE  10  m 

central  Panama 
156  million  cells 

SRTM  30  m 
all  of  Panama 

241  million  cells 

old  r.watershed , 
all  in  memory 

23.2min 

23  h  1 0  min 

30  h  46  min 

>  2  days 
not  measured 

>  2  days 
not  measured 

new  r.watershed , 
all  in  memory 

0.07  min 

0.83  min 

1.05  min 

9.97  min 

12.3  min 

new  r.watershed , 
data  on  disk 

0.1 5  min 

1.73  min 

2.05  min 

30.7  min 

32.18  min 

r.terraflow 

0.45  min 

4.4  min 

5.05  min 

66.95  min 

58.7  min 

Table  2.  Median  distances  of  reference  points  to  nearest  extracted 
stream  in  meters. 


GPS 

filling 

IRA 

LCP 

Landsat 

filling 

IRA 

LCP 

IFSARE  10  m 

62.65 

N/A* 

34.40 

43.75 

N/A* 

39.25 

IFSARE  30  m 

77.15 

71.42 

70.71 

48.88 

48.27 

47.75 

SRTM  30  m 

99.78 

92.46 

89.13 

84.25 

82.81 

83.72 

SRTM  90  m 

124.84 

119.69 

121.00 

94.30 

94.69 

94.92 

*  Computation  cancelled  by  software  with  out  of  memory  error. 


IFSARE  IFSARE  SRTM  SRTM 

10m  30m  30m  90m 


□  GPS  points 
Q]  Landsat  points 


10m  30m  30m  90m 


I  I  GPS  points 
|  |  Landsat  points 


Fig.  5.  Mean  and  standard  deviation  of  the  distance  differences  be¬ 
tween  sink  filling  and  least-cost  path  search.  A  positive  distance  dif¬ 
ference  indicates  that  streams  extracted  with  least-cost  path  search 
are  closer  to  reference  points  than  streams  extracted  with  sink  fill¬ 
ing. 


Fig.  6.  Mean  and  standard  deviation  of  the  distance  differences 
between  impact  reduction  approach  and  least-cost  path.  A  posi¬ 
tive  distance  difference  indicates  that  streams  extracted  with  least- 
cost  path  search  are  closer  to  reference  points  than  streams  extracted 
with  the  impact  reduction  approach. 


lize,  hence  the  out-of-memory  error  was  regarded  as  a  lim¬ 
itation  of  this  application  and  not  as  a  hardware  limitation. 
Therefore,  TAS  was  used  only  for  the  sink  treatment  meth¬ 
ods  in  comparison  with  the  30  m  and  90  m  DEM  covering 
the  central  Panama  subregion  (Fig.  2)  and  results  of  the  IRA 
method  were  not  available  for  10  m  IFSARE.  At  the  time  of 
writing,  IRA  was  not  included  in  WhiteboxGAT  (Lindsay, 
2009)  that  has  replaced  TAS. 

3.2  Comparison  between  different  methods  of 
sink  treatment 


Panama  exceeded  available  system  memory.  TAS  GIS  is  a 
32  bit  application  for  Microsoft  Windows  and  can  only  uti¬ 
lize  as  much  memory  as  a  32  bit  application  can  manage, 
i.e.  3  GB.  As  opposed  to  r.terraflow  and  r.watershed ,  TAS 
GIS  keeps  all  intermediate  data  in  memory,  thus  limiting  the 
size  of  the  DEMs  that  can  be  processed.  The  test  system  pro¬ 
vided  as  much  memory  as  TAS  GIS  could  theoretically  uti- 


To  evaluate  the  impact  of  different  approaches  in  the  han¬ 
dling  of  elevation  surface  depressions,  median  distances  of 
reference  points  to  extracted  drainage  paths  for  each  sink 
treatment  method  and  each  processed  DEM  were  computed 
and  the  results  are  summarized  in  Table  2.  As  expected,  the 
location  of  extracted  drainage  paths  improved  with  increas¬ 
ing  detail  available  in  the  DEM.  While  the  drainage  paths  for 
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GPS  field  points  Landsat 


□  lcp 

1  |  Sink  filling 

□  IRA 


Fig.  7.  Mean  and  standard  deviation  of  the  distance  differences 
between  the  30  m  SRTM  DEM  and  the  30  m  IFSARE  DEM.  A  pos¬ 
itive  distance  difference  indicates  that  streams  extracted  from  the 
30  m  IFSARE  DEM  are  closer  to  reference  points  than  streams  ex¬ 
tracted  from  the  30  m  SRTM  DEM.  LCP:  least-cost  path  search; 
IRA:  impact  reduction  approach. 


GPS  field  points  Landsat 


□  lcp 

I  |  Sink  filling 

□  IRA 


Fig.  8.  Mean  and  standard  deviation  of  the  distance  differences  be¬ 
tween  the  90  m  and  the  30  m  SRTM  DEM.  A  positive  distance  dif¬ 
ference  indicates  that  streams  extracted  from  the  30  m  SRTM  DEM 
are  closer  to  reference  points  than  streams  extracted  from  the  90  m 
SRTM  DEM.  LCP:  least-cost  path  search;  IRA:  impact  reduction 
approach. 


the  90  m  SRTM  DEM  were  about  100  m  away  from  refer¬ 
ence  points,  the  median  distances  improved  down  to  34.4  m 
for  the  IFSARE  10  m  DEM. 

Drainage  paths  extracted  from  the  LCP  flow  accumulation 
were  closer  to  the  GPS  field  points  (t- tests,  all  t^i  <  —  5.4,  all 
p  <  0.0001)  for  all  DEMs  at  all  tested  resolutions  (IFSARE 
10  m,  30  m,  and  SRTM  30  m,  90  m)  when  compared  with  the 
sink  filling  method  (Fig.  5  showing  bar  graphs  with  mean  and 
standard  deviation  of  the  distance  differences).  The  drainage 
paths  extracted  from  LCP  were  also  closer  to  the  points  dig¬ 
itized  from  Landsat  for  all  DEMs  (t- tests,  all  £994  <  —2.3, 
all  p  <  0.0001),  except  for  the  90  m  SRTM  where  the  dis¬ 
tance  difference  was  not  significant  for  (t-test,  £994  =  —  1 .213, 
p  =  0.23). 

When  compared  to  the  IRA  method  (Fig.  6  showing  bar 
graphs  with  mean  and  standard  deviation  of  the  distance  dif¬ 
ferences),  drainage  paths  extracted  from  LCP  flow  accumula¬ 
tion  were  closer  to  the  GPS  field  points  for  the  30  m  IFSARE 


(t-test,  ^337  =  —2.361,  p  =  0.018)  and  30  m  SRTM  (t-test, 
£337  =  —3.088,  p  =  0.002),  but  for  the  90m  SRTM  the  dif¬ 
ference  was  not  significant  (t-test,  £337  =  —1.745,  p  =  0.08). 
For  the  points  digitized  from  Landsat,  drainage  paths  ex¬ 
tracted  from  LCP  accumulation  were  not  closer  than  IRA 
drainage  paths  for  any  of  the  tested  DEMs  (t-tests,  all  ab¬ 
solute  £994  <  1.0,  all  p  >  0.3). 

3.3  Comparison  between  different  levels  of  detail 
at  constant  resolution 

The  10  m  IFSARE  and  90  m  SRTM  were  interpolated  to  the 
same  resolution  of  30  m.  We  investigated  how  much  the  dif¬ 
ferent  levels  of  detail  at  identical  resolution  influenced  the  ac¬ 
curacy  of  extracted  drainage  paths  for  all  three  sink  treatment 
methods  (Fig.  7  showing  bar  graphs  with  mean  and  standard 
deviation  of  the  distance  differences).  The  distance  between 
the  reference  points  and  drainage  paths  extracted  from  the 
30  m  IFSARE  DEM  was  subtracted  from  the  distance  be¬ 
tween  the  reference  points  and  drainage  paths  extracted  from 
the  30  m  SRTM  DEM  to  measure  the  difference  in  accuracy 
of  drainage  paths  extracted  from  these  two  DEMs.  For  larger 
rivers  in  flat  areas  (Landsat  points),  the  accuracy  of  drainage 
path  locations  was  considerably  higher  for  the  30  m  IFSARE 
DEM  than  for  the  30  m  SRTM  DEM  (t-test,  all  ^994  <  — 13.0, 
all  p  <  0.001;  on  average  35  m  closer)  for  all  methods.  For 
smaller  rivers  not  broader  than  30  m  in  mountainous  regions 
(GPS  points),  only  IRA  showed  a  significant  improvement 
in  accuracy  from  the  30  m  SRTM  DEM  to  the  30  m  IFSARE 
DEM  (t-test,  £337  =  —3.315,  p  =  0.001,  all  other  t^i  >  — 1.9 
and  <0,  p  >  0.05).  However,  as  shown  in  the  previous  sec¬ 
tion  for  the  30  m  IFSARE  DEM,  IRA  drainage  path  locations 
in  mountainous  regions  were  significantly  less  accurate  than 
LCP  drainage  path  locations. 

3.4  Effects  of  RST  resampling  on  the  accuracy  of 
drainage  network  extraction 

In  order  to  assess  the  effect  of  SRTM  resampling  from  90  m 
to  30  m  resolution  by  the  RST  method,  the  drainage  paths 
extracted  from  the  90  m  SRTM  DEM  were  compared  with 
drainage  paths  extracted  from  the  resampled  30  m  SRTM 
DEM  by  calculating  the  difference  in  drainage  path  distances 
to  reference  points  (Fig.  8  showing  bar  graphs  with  mean  and 
standard  deviation  of  the  distance  differences).  Correspond¬ 
ing  median  distances  of  reference  points  to  drainage  path  lo¬ 
cations  are  given  in  Table  2.  The  accuracy  of  drainage  path 
locations  was  improved  by  resampling  for  larger  rivers  in  flat 
landscape  for  all  sink  treatment  methods  (t-tests,  all  £994  < 
—6.5,  all  p  <  0.001),  however,  for  GPS  points  in  mountain¬ 
ous  regions  only  the  results  obtained  by  the  least-cost  path 
search  were  improved  (t-test,  £337  =  —3.734,  p  <  0.001,  all 
other  £337  >  —  1  and  <0,  p  >  0.3). 
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Fig.  9.  A  typical  example  for  the  impact  of  sink  filling  along  a  section  of  a  river  course  extracted  from  the  30  m  SRTM  DEM.  Along  this 
section,  92%  of  all  elevation  values  were  modified.  Parts  of  the  river  course  obscured  by  trees  are  highlighted  in  the  GoogleEarth  sreenshot. 


Fig.  10.  Typical  effects  of  different  sink  treatment  methods  on  extracted  stream  courses  in  (A)  mountainous  areas  and  (B)  costal  plains.  Filled 
sinks  are  grayed.  Streams  extracted  from  sink-filling  in  red  are  overlaid  by  streams  extracted  from  least-cost  path  search  in  blue.  Points  of 
confluences  (yellow  cross)  in  mountainous  areas  can  be  shifted  considerably,  and  in  coastal  plains,  sink  filling  can  lead  to  long  straight  lines. 
A  shaded  relief  of  the  10  m  IFSARE  DEM  was  used  as  background. 


4  Discussion 

Digital  Terrain  Models,  as  an  approximation  of  real  terrain  at 
a  given  level  of  detail,  include  random  noise  and  errors.  For 
example,  an  in  depth  analysis  of  SRTM  accuracy  and  source 
of  errors  is  provided  by  Rodriguez  et  al.  (2006)  and  Farr  et 
al.  (2007).  The  errors  can  have  significant  impact  on  the  re¬ 
sults  of  DEM  analysis  and  their  treatment  has  been  the  focus 
of  broad  research  (e.g.  Jenson  and  Domingue,  1988;  Lindsay 
and  Creed,  2005;  Hutchinson,  2000).  Sink  (depression)  re¬ 
moval  was  designed  to  remove  small  depressions  introduced 
as  artifacts  of  elevation  data  processing  for  the  first  genera¬ 
tion  of  DEMs  to  facilitate  continuous  flow  routing  (Jenson 
and  Domingue,  1988).  This  procedure,  referred  to  as  hydro- 
logical  conditioning,  has  become  widespread  even  for  next 
generations  of  elevation  data  and  is  now  a  pre-processing 


step  required  by  a  number  of  applications  for  hydrological 
analysis  and  modeling.  In  order  to  alleviate  the  problem  of 
routing  surface  flow  through  the  flat  areas  created  by  sink 
filling,  efforts  have  been  made  to  add  a  gradient  to  these  flat 
areas  (e.g.  Martz  and  Garbrecht,  1998;  Wang  and  Liu,  2006; 
Grimaldi  et  al.,  2007;  Santini  et  al.,  2009)  Nevertheless,  as 
modem  high  resolution  data  such  as  LiDAR  demonstrate, 
actual  terrain  includes  many  tme  depressions,  especially  in 
riffle-and-pool  portions  of  mountainous  rivers,  and  on  flood- 
plains  and  coastal  plains  (Notebaert  et  al.,  2009).  There¬ 
fore,  sink  filling  is  not  always  an  appropriate  approach  for 
hydrological  conditioning  of  a  DEM  because  the  measured 
elevation  data  are  replaced  with  new  elevation  values  based 
on  the  rather  unrealistic  assumption  that  the  terrain  has  no 
depressions  (Jenson  and  Domingue,  1988;  Tarboton,  1997), 
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be  they  artificial  or  real.  Although  many  hydrologic  mod¬ 
els  require  depression-less  DEMs,  removal  of  sinks  can  alter 
the  true  hydrological  state  that  includes  standing  water  in  de¬ 
pressions.  Moreover,  in  IFSARE  and  SRTM  DEMs  that  in¬ 
clude  vegetation  cover,  many  of  the  observed  depressions  are 
caused  by  gaps  in  canopy,  canopy  closing  over  river  courses, 
or  by  vegetation  partially  obstructing  narrow  and  deep  val¬ 
leys.  The  filling  of  such  features  requires  altering  elevations 
by  several  meters  and  over  areas  of  hundreds  of  m2,  leading 
to  large  flats  and  an  artificial  drainage  network  geometry  that 
does  not  represent  reality.  Figure  9  presents  an  example  of 
the  barriers  across  the  river  channel  created  by  overhanging 
trees  and  a  profile  that  illustrates  the  extensive  modification 
in  elevation  that  can  result  from  such  filling.  Preserving  most 
of  the  drainage  directions  of  the  original  DEM  is  preferable 
to  hydrological  conditioning  in  this  case. 

Our  results  suggest  that  the  adverse  effects  of  sink  filling 
become  more  pronounced  for  higher  resolution  DEMs  (see 
Fig.  10  for  drainage  paths  extracted  from  10  m  IFSARE). 
By  contrast,  at  lower  resolutions  such  as  90  m,  the  coarsest 
resolution  tested  in  this  study,  there  was  often  no  significant 
difference  in  drainage  path  locations  between  the  compared 
methods,  particularly  for  Landsat  reference  points.  Land- 
sat  points  where  predominantly  located  in  flood  plain  and 
coastal  plain  landscapes  where  both  the  horizontal  and  the 
vertical  resolution  of  90  m  SRTM  were  often  not  sufficient  to 
locate  rivers  and  no  method  was  able  to  accurately  delineate 
drainage  paths  in  these  areas  from  90  m  SRTM  (Hancock  et 
al.,  2006;  Valeriano  et  ah,  2006). 

When  comparing  SRTM  and  IFSARE  DEMs  at  the  same 
resolution  of  30  m,  the  drainage  networks  extracted  from 
the  IFSARE  DEM  were,  as  expected,  more  accurate  than 
drainage  networks  extracted  from  SRTM  DEM,  particularly 
for  Landsat  points  in  relatively  flat  areas.  As  shown  in  previ¬ 
ous  studies  (Hancock  et  al.,  2006;  Valeriano  et  ah,  2006), 
the  SRTM  DEM  close  to  its  native  resolution  of  3  arc  sec 
does  not  provide  sufficient  vertical  detail  to  determine  river 
courses  in  flat  areas.  The  applied  RST  resampling  method 
significantly  improved  the  accuracy  of  extracted  drainage 
path  locations  for  SRTM  in  regions  with  low  topography. 
Apparently  the  level  of  detail  has  been  not  just  increased  but 
also  improved  because  resampling  to  a  higher  horizontal  (and 
for  SRTM  also  higher  vertical)  resolution  seems  to  be  help¬ 
ful  in  improving  the  results  of  watershed  analysis  from  low- 
resolution  DEMs  (see  also  Valeriano  et  ah,  2006).  The  global 
SRTM  coverage  at  3  arc  sec  horizontal  resolution  has  been 
created  by  averaging  the  original  1  arc  sec  coverage  (Farr  et 
al.,  2007).  Depending  on  the  method  used,  reinterpolating 
the  3  arc  sec  version  back  to  1  arc  sec  can  result  in  narrowed 
valleys  and  ridges  and  increased  channel  sinuosity.  Interest¬ 
ingly,  only  the  LCP  search  method  benefited  from  the  reinter¬ 
polation  in  mountainous  regions  and  provided  here  the  most 
accurate  drainage  path  locations  for  the  30  m  SRTM  DEM. 

It  is  important  to  note  that  the  differences  in  the  tested 
methods  were  significant  due  to  the  particularly  challenging 


study  region.  Smaller  differences  between  the  methods  could 
be  expected  in  topography  with  gentle  slopes  and  uniformly 
low  vegetation  cover.  Although  network  connectivity,  topol¬ 
ogy  and  channel  morphometry  has  not  been  investigated,  the 
visual  inspection  of  the  results  indicates  that  there  were  dif¬ 
ferences  between  the  results  obtained  by  different  methods 
and  resolutions  in  the  coastal  plain  areas.  None  of  the  tested 
DEMs  or  methods  produced  sufficiently  reliable  results  and 
higher  resolution  did  not  always  lead  to  correct  topology. 
These  results  are  similar  to  those  found  for  LiDAR-based 
DEMs  (Colson,  2006)  and  require  further  research. 

Although  the  tests  presented  here  demonstrated  the  per¬ 
formance  of  the  LCP  method  for  radar-based  data,  the  new, 
improved  version  of  r.  watershed  has  in  fact  been  developed 
using  LiDAR  data  provided  by  the  USGS  in  addition  to  the 
1  m  LiDAR  DEM  examples  and  data  provided  by  Neteler  and 
Mitasova  (2008).  The  presented  method  processes  a  DEM 
without  modifying  its  elevation  values  and  is  thus  to  a  de¬ 
gree  (quantitatively  narrowed  by  this  study)  dependent  on 
the  reliability  of  the  provided  DEM.  The  LCP  search  ap¬ 
proach,  despite  its  capability  to  route  flow  through  depres¬ 
sions,  will  therefore,  like  other  methods,  fail  to  accurately 
extract  drainage  networks  if  the  input  DEM  deviates  too  far 
from  real  topography. 

5  Conclusions 

We  presented  a  method  for  fast  computation  of  flow  accu¬ 
mulation  and  drainage  network  extraction  for  large  DEMs 
with  nested  depressions  and  evaluated  it  against  other  com¬ 
monly  used  methods  for  sink  treatment.  The  accuracy  assess¬ 
ment  using  ground  control  points  and  Landsat  satellite  im¬ 
agery  provided  insight  into  the  accuracy  of  drainage  network 
extraction  from  radar  elevation  data  in  a  challenging  tropi¬ 
cal  forest  environment  and  coastal  plain  setting.  The  results 
suggest  that  the  conceptually  simple  sink  filling  approach  is 
not  suitable  for  the  IFSARE  and  SRTM  DEMs  in  regions 
with  significant  vegetation  cover  whereas  both  the  tested  im¬ 
pact  reduction  approach  of  Lindsay  and  Creed  (2005)  and  the 
LCP  search  provide  more  accurate  drainage  networks. 

The  performance  testing  has  demonstrated  that  the  new 
implementation  dramatically  improves  computational  effi¬ 
ciency  while  preserving  the  high  accuracy  of  the  LCP  rout¬ 
ing  capabilities.  The  increase  in  computational  performance 
is  particularly  relevant  for  the  modem  mapping  technologies 
that  can  rapidly  produce  massive  DEMs  at  high  resolutions 
for  large  areas  and  the  slow  processing  and  analysis  has  be¬ 
come  bottleneck  in  their  application.  Fast  processing  of  mas¬ 
sive  DEMs  opens  their  application  to  a  new  level  of  detail  and 
spatial  extent,  for  example  in  rapid  response  operations  or  for 
mapping  in  remote  regions  where  the  streams  are  covered  by 
dense  vegetation  and  no  reliable  stream  data  are  available. 


50 


Acknowledgements.  The  support  of  the  US  Army  Research  Office 

for  this  research  is  gratefully  acknowledged.  We  would  also  like  to 

thank  two  anonymous  reviewers  for  their  valuable  comments. 

Edited  by:  R.  Purves 

References 

Andersen,  H.  E.,  Reutebuch,  S.  E.,  and McGaughey,  R.  J.:  Accuracy 
of  an  IFSAR-derived  digital  terrain  model  under  a  conifer  forest 
canopy,  Can.  J.  Remote  Sens.,  31,  283-288,  2005. 

Arge,  L.,  Chase,  J.  S.,  Halpin,  P.  N.,  Toma,  L.,  Vitter,  J.  S.,  Urban, 
D.,  and  Wickremesinghe,  R.:  Flow  computation  on  massive  grid 
terrains,  Geolnformatica,  7,  283-313,  2003. 

Atkinson,  M.  D.,  Sack,  J.-R.,  Santoro,  N.,  and  Strothotte,  T.:  Min- 
max  heaps  and  generalized  priority  queues,  Programming  tech¬ 
niques  and  Data  structures,  Comm.  ACM,  29,  996-1000,  1986. 

Colson,  T.  P:  Stream  network  delineation  from  high-resolution 
digital  elevation  models,  Ph.D.  Dissertation,  Department  of 
Forestry  &  Environmental  Resources,  North  Carolina  State  Uni¬ 
versity,  Raleigh,  NC,  available  at  http://www.lib.ncsu.edu/theses/ 
available/etd- 10302006- 122024/,  2006. 

Costa-Cabral,  M.  C.  and  Burges,  S.  J.:  Digital  elevation  model  net¬ 
works  (DEMON):  A  model  of  flow  over  hillslopes  for  computa¬ 
tion  of  contributing  and  dispersal  areas,  Water  Resour.  Res.,  30, 
1681-1692,  1994. 

Danner,  A.,  Yi,  K.,  Moelhave,  T.,  Agarwal,  P.  K.,  Arge,  L.,  and  Mi- 
tasova,  FL:  TerraStream:  From  Elevation  Data  to  Watershed  Hi¬ 
erarchies,  Proc.  ACM  GIS,  28,  doi:  10. 1 145/1341012. 1341049, 
2007. 

Ehlschlaeger,  C.:  Using  the  A*  Search  Algorithm  to  Develop  Hy¬ 
drologic  Models  from  Digital  Elevation  Data,  Proceedings  of  In¬ 
ternational  Geographic  Information  Systems  (IGIS)  Symposium, 
Baltimore,  MD,  USA,  275-281,  1989. 

Fairfield,  J.  and  Leymarie,  P:  Drainage  networks  from  grid  digital 
elevation  models,  Water  Resour.  Res.,  27,  709-717,  1991. 

Farr,  T.  G.,  Rosen,  A.  R.,  Caro,  E.,  Crippen,  R.,  Duren,  R., 
Hensley,  S.,  Kobrick,  M.,  Paller,  M.,  Rodriguez,  E.,  Roth, 

L. ,  Seal,  D.,  Shaffer,  S.,  Shimanda,  J.,  Umland,  J.,  Werner, 

M. ,  Oskin,  M.,  Burbank,  D.,  and  Alsdorf,  D.:  The  Shut¬ 
tle  Radar  Topography  Mission,  Rev.  Geophys.,  45,  RG2004, 
doi:  10.1 029/2005RG000 1 83 ,  2007. 

Freeman,  T.  G.:  Calculating  catchment  area  with  divergent  flow 
based  on  a  regular  grid,  Comput  Geosci.,  17,  413^122,  1991. 

Garbrecht,  J.  and  Martz,  L.  W.:  The  assignment  of  drainage  direc¬ 
tion  over  flat  surfaces  in  raster  digital  elevation  models,  J.  Hy- 
drol.,  192,204-213,  1997. 

GRASS  Development  Team:  Geographic  Resources  Analysis  Sup¬ 
port  System  (GRASS)  Software,  Version  6.3.0.  http://grass. 
osgeo.org/download/software_old.php#g630,  2008. 

GRASS  Development  Team:  Geographic  Resources  Analysis  Sup¬ 
port  System  (GRASS)  Software,  Version  6.4.0,  Open  Source 
Geospatial  Foundation,  http://grass.osgeo.org,  2010. 

Grimaldi,  S.,  Nardi,  F.,  Di  Benedotto,  F.,  Istanbulluoglu,  E.,  and 
Bras,  R.  L.:  A  physically-based  method  for  removing  pits  in  dig¬ 
ital  elevation  models,  Adv.  Water  Resour.,  30,  2151-2158,  2007. 

Hancock,  G.  R.,  Martinez,  C.,  Evans,  K.  G.,  and  Moliere,  D.  R.:  A 
comparison  of  SRTM  and  high-resolution  digital  elevation  mod¬ 
els  and  their  use  in  catchment  geomorphology  and  hydrology: 


Australian  examples,  Earth  Surf.  Proc.  Land.,  31,  1394-1412, 
2006. 

Hart,  P.  E.,  Nilsson,  N.  J.,  and  Raphael,  B.:  A  Formal  Basis  for  the 
Heuristic  Determination  of  Minimum  Cost  Paths,  IEEE  T.  Syst. 
Sci.  Cyb.,4,  100-107,  1968. 

Holmgren,  P.:  Multiple  flow  direction  algorithms  for  runoff  mod¬ 
eling  in  grid  based  elevation  models:  An  empirical  evaluation, 
Hydrol.  Process.,  8,  327-334,  1994. 

Jenson,  S.  K.  and  Domingue,  J.  O.:  Extracting  topographic  struc¬ 
ture  from  digital  elevation  data  for  geographic  information  sys¬ 
tem  analysis,  Photogramm.  Eng.  Rem.  S.,  54(11),  1593-1600, 
1988. 

Kinner,  D.,  Mitasova,  H.,  Harmon,  R.  S.,  Toma,  L.,  and  Stal- 
lard,  R.:  GIS-based  Stream  Network  Analysis  for  The  Chagres 
River  Basin,  Republic  of  Panama.  In:  Harmon  R  (ed)  The  Rio 
Chagres:  A  Multidisciplinary  Profile  of  a  Tropical  Watershed, 
Springer/Kluwer,  83-95,  2005. 

Lehner,  B.  and  Doll,  P.:  Development  and  validation  of  a  global 
database  of  lakes,  reservoirs  and  wetlands,  J.  Hydrol.,  296,  1-22, 
2004. 

Lindsay,  J.  B.:  Whitebox  http://www.uoguelph.ca/~hydrogeo/ 
Whitebox/index.html,  2009. 

Lindsay,  J.  B.  and  Creed,  F.:  Removal  of  artefact  depressions  from 
digital  elevation  models:  towards  a  minimum  impact  approach, 
Hydrol.  Process.,  19,  3113-3126,  2005. 

Martz,  L.  W.  and  Garbrecht,  J.:  The  treatment  of  flat  areas  and  de¬ 
pressions  in  automated  drainage  analysis  of  raster  digital  eleva¬ 
tion  models,  Hydrol.  Process.,  12,  843-855,  1998. 

Metz,  M.  and  Ehlschlaeger,  C.:  Watershed  analysis  program 
r.watershed,  source  code,  https://trac.osgeo.org/grass/browser/ 
grass/branches/releasebranch_6_4/raster/r. watershed,  20 1 0. 

Mitasova,  H.  and  Mitas,  L.:  Interpolation  by  Regularized  Spline 
with  Tension:  I.  Theory  and  Implementation,  Math.  Geol.,  25, 
641-655,  1993. 

Neteler,  M.  and  Mitasova,  H.:  Open  Source  GIS:  A  GRASS  GIS 
Approach,  Third  Edition,  Springer  New  York,  406  pp.,  2008. 

North  Carolina  Division  of  Water  Quality  (NCWQ):  Methodol¬ 
ogy  for  Identification  of  Intermittent  and  Perennial  Streams 
and  Their  Origins,  Version  4.1.1,  http://portal.ncdenr.org/web/ 
wq/swp/ws/401/waterresources/streamdeterminations,  Effective 
Date:  1  September  2010. 

Notebaert,  B.,  Verstraeten,  G.,  Govers,  G.,  and  Poesen,  J:  Qual¬ 
itative  and  quantitative  applications  of  LiDAR  imagery  in  flu¬ 
vial  geomorphology,  Earth  Surf.  Proc.  Landforms,  34,  217-231, 
2009. 

O’Callaghan,  J.  F.  and  Mark,  D.  M.:  The  Extraction  of  Drainage 
Networks  from  Digital  Elevation  Data,  Computer  Vision,  Graph¬ 
ics,  and  Image  Processing,  28,  223-344,  1984. 

Peckham,  S.  D.:  Efficient  extraction  of  river  networks  and  hydro- 
logic  measurements  from  digital  elevation  data,  in  Barndorff- 
Nielsen  and  others,  eds.,  Stochastic  Methods  in  Hydrology: 
Rain,  Landforms  and  Floods:  Singapore,  World  Scientiifc,  173— 
203,  1998. 

Planchon,  O.  and  Darboux,  F.:  A  fast,  simple  and  versatile  algo¬ 
rithm  to  fill  the  depressions  of  digital  elevation  models,  Catena, 
46,  159-176,2001. 

Quinn,  R,  Beven,  K.,  Chevallier,  P.,  and  Planchon,  O.:  The  predic¬ 
tion  of  hillslope  flow  paths  for  distributed  hydrological  modelling 
using  digital  terrain  models,  Hydrol.  Process.,  5,  59-79,  1991. 


51 


Quinn,  P.  R,  Beven,  K.  J.,  and  Lamb,  R.:  The  ln(a/tan/B)  index: 
How  to  calculate  it  and  how  to  use  it  within  the  topmodel  frame¬ 
work,  Hydro!  Process.,  9,  161-182,  1995. 

Rieger,  W.:  A  phenomenon-based  approach  to  upslope  area  and  de¬ 
pressions  inDEMs,  Hydro!  Process.,  12,  857-872,  1998. 

Rivix  Limited  Liability  Company:  RiverTools™User’s  Guide  re¬ 
lease  2001,  Boulder,  CO,  Research  Systems,  Inc.,  202  pp.,  2001. 

Rodriguez,  E.,  Morris,  C.  S.,  and  Belz,  J.  E.:  A  global  assessment 
of  the  SRTM  performance,  Photogramm.  Eng.  Rem.  S.,  72,  249- 
260,  2006. 

Santini,  M.,  Grimaldi,  S.,  Rulli,  M.  C.,  Petroselli,  A.,  and  Nardi,  F.: 
Pre-processing  algorithms  and  landslide  modeling  on  remotely 
sensed  DEMs,  Geomorphology,  113,  110-125,  2009. 


Tarboton,  D.  G.:  A  New  Method  for  the  Determination  of  Flow  Di¬ 
rections  and  Contributing  Areas  in  Grid  Digital  Elevation  Mod¬ 
els,  Water  Resour.  Res.,  33,  309-319,  1997. 

Valeriano,  M.  M.,  Kuplich,  T.  M.,  Storino,  M.,  Amaral,  B.  D., 
Mendes,  J.  N.,  and  Lima,  D.  J.:  Modeling  small  watersheds  in 
Brazilian  Amazonia  with  shuttle  radar  topographic  mission  90m 
data,  Comput.  Geosci.,  32,  1169-1181,  2006. 

Wang,  L.  and  Liu,  H.  An  efficient  method  for  identifying  and  fill¬ 
ing  surface  depressions  in  digital  elevation  models  for  hydrologic 
analysis  and  modelling,  Int.  J.  Geogr.  Inf.  Sc!,  20,  193-213, 
2006. 

World  Wildlife  Fund,  HydroSHEDS,  http://hydrosheds.cr.usgs.gov/ 
index.php,  2009. 


